{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "19095b39-daba-4632-8b23-35b22ab71a8c",
   "metadata": {
    "tags": []
   },
   "source": [
    "# 快速傅里叶变换fft\n",
    "具体要求参考[链接](https://www.wolai.com/zzutai/ne85AbnL15kmZPPJRQ2q4j)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "4e42cc3e-b29d-463d-aa51-5a1dae08a51b",
   "metadata": {
    "tags": []
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAD4CAYAAACg9uHUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAesklEQVR4nO3df5DkdX3n8ec7y5JbkbvVMBAYSMBzmQqJFzehIBSV1JDkBLZMWL3oQVIRL9at5mJVsJQEohVNUndyoczlUhjIqhR4MQiVwErixpUT5tASRGAXFsQNCxrc2T22UAcZmdNled8f853YO3Rvz/b32/P9fnuej6qu7v7+6O+rZz7f97y/3f3ticxEkiRJ9fihugNIkiStZDZjkiRJNbIZkyRJqpHNmCRJUo1sxiRJkmp0VN0Bulm7dm2++tWvrjtGX9/97nc55phj6o7Rlzmr15asbcn5wAMPPJOZY3XnqIL1q1ptyQntyWrO6pWtYY1sxk444QTuv//+umP0NTU1xeTkZN0x+jJn9dqStS05I+Kf685QFetXtdqSE9qT1ZzVK1vDfJtSkiSpRjZjkiRJNbIZkyRJqpHNmCRJUo1sxiRJkmrUtxmLiFMi4q6IeCwiHo2I3y2mvzIi7oiIx4vrV/RY/4KI2BURuyPiiqqfwJHasn2ac6+6k9Ou+DTnXnUnW7ZP1x1JaqRR2FdGrX5JGk1LeWXsBeDdmfkTwM8BvxMRZwBXAJ/LzHXA54r7h4iIVcCHgQuBM4BLinVrsWX7NFfeupPpmTkSmJ6Z48pbd7byj4w0TCO0r4xM/ZI0uvo2Y5m5LzMfLG4/BzwGjAMXATcWi90IbOyy+lnA7sx8MjO/D3yyWK8WV2/bxdyBg4dMmztwkKu37aopkdRMo7KvjFL9kjS6juhLXyPiVGA98CXghMzcB/MFLyKO77LKOPCNjvt7gLN7PPYmYBPA2NgYU1NTRxJtSaZn5npOH2R7s7OzQ8lZNXNWry1ZB81Z9b7SBG2vX1Ub9TFch7ZkNWfzLLkZi4iXA38HXJaZ34mIJa3WZVp2WzAzNwObASYmJnIY37o7fu+dXf/IjK9dM9C3/Lbl24HNWb22ZB00Z9X7St1GoX5VbdTHcB3aktWczbOksykjYjXzhewTmXlrMfnpiDixmH8isL/LqnuAUzrunwzsHTxuOZefP8Ga1asOmbZm9SouP3+ipkRSM43SvjIq9UvS6FrK2ZQBfAx4LDP/rGPW7cClxe1LgU91Wf3LwLqIOC0ijgYuLtarxcb143zwja/h6FXzT3t87Ro++MbXsHH9eF2RpEYalX1llOqXpNG1lLcpzwV+E9gZETuKaX8AXAXcEhFvA54C3gQQEScBH83MDZn5QkS8E9gGrAKuz8xHK34OR2Tj+nFuuu8pAG5++zl1RpEabUT2lZGqX5JGU99mLDO/QPfPTgD8Upfl9wIbOu5vBbYOGlCSBmX9ktQGfgO/JElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1OqrfAhFxPfB6YH9m/lQx7WZgolhkLTCTma/tsu7XgeeAg8ALmXlmJaklSZJGRN9mDLgBuAb4+MKEzPyPC7cj4kPAs4dZ/7zMfGbQgJJUhgeUkpqubzOWmXdHxKnd5kVEAG8GfrHiXJJUlRvwgFJSg5X9zNjPA09n5uM95ifw2Yh4ICI2ldyWJB2xzLwb+Fa3eR0HlDctayhJ6rCUtykP5xIOX8TOzcy9EXE8cEdEfLUojC9RNGubAMbGxpiamioZrbeZmTmA0tuYnZ0das6qmLN6bclaNmdV+0qDLfWAMoG/yszN3RZazvpVlZUyhpdTW7Kas4Eys+8FOBV4ZNG0o4CngZOX+BgfAN6zlGVPP/30HKY3X/fFfPN1Xyz9OHfddVf5MMvAnNVrS9ayOavaV/oB7s8l1IZBL91qWDH9WuDdh1nvpOL6eOAh4Bf6bWvY9asqK2UML6e2ZDVn9crWsDJvU/4y8NXM3NNtZkQcExHHLtwGXgc8UmJ7klSZiDgKeCNwc69lMnNvcb0fuA04a3nSSVpJ+jZjEXETcA8wERF7IuJtxayLWfQWZUScFBFbi7snAF+IiIeA+4BPZ+ZnqosuSaV4QCmpEZZyNuUlPaa/tcu0vcCG4vaTwE+XzCdJpRQHlJPAcRGxB3h/Zn6MHgeUwEczcwPzB5S3zX/Gn6OAv/GAUtIwlP0AvyQ1mgeUkprOf4ckSZJUI5sxSZKkGtmMSZIk1chmTJIkqUY2Y5IkSTWyGZMkSaqRzZgkSVKNbMYkSZJqZDMmSZJUI5sxSZKkGtmMSZIk1chmTJIkqUY2Y5IkSTWyGZMkSaqRzZgkSVKN+jZjEXF9ROyPiEc6pn0gIqYjYkdx2dBj3QsiYldE7I6IK6oMLkmSNAqW8srYDcAFXab/j8x8bXHZunhmRKwCPgxcCJwBXBIRZ5QJK0lHygNKSU3XtxnLzLuBbw3w2GcBuzPzycz8PvBJ4KIBHkeSyrgBDyglNdhRJdZ9Z0S8BbgfeHdmfnvR/HHgGx339wBn93qwiNgEbAIYGxtjamqqRLTDm5mZAyi9jdnZ2aHmrIo5q9eWrGVzVrWv1Ckz746IUwdY9V8OKAEiYuGA8isVxpOkgZuxa4E/AbK4/hDwW4uWiS7rZa8HzMzNwGaAiYmJnJycHDBaf9fuugeAyclzSj3O1NQUw8xZFXNWry1Zy+asal9pqMoOKJfzYLIqK+WAYjm1Jas5m2egZiwzn164HREfAf6hy2J7gFM67p8M7B1ke5JUsUoPKJfzYLIqK+WAYjm1Jas5m2egr7aIiBM77r4BeKTLYl8G1kXEaRFxNHAxcPsg25OkKmXm05l5MDNfBD7C/FuSi3lAKWlZLOWrLW4C7gEmImJPRLwN+NOI2BkRDwPnAe8qlj0pIrYCZOYLwDuBbcBjwC2Z+eiQnockLZkHlJKapO/blJl5SZfJH+ux7F5gQ8f9rcBLzlKSpOVSHFBOAsdFxB7g/cBkRLyW+bcdvw68vVj2JOCjmbkhM1+IiIUDylXA9R5QShqGMmdTSlLjeUApqen8d0iSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNWobzMWEddHxP6IeKRj2tUR8dWIeDgibouItT3W/XpE7IyIHRFxf4W5JWlJrGGSmm4pr4zdAFywaNodwE9l5r8D/gm48jDrn5eZr83MMweLKEml3IA1TFKD9W3GMvNu4FuLpn02M18o7t4LnDyEbJJUmjVMUtNFZvZfKOJU4B8y86e6zPt74ObM/Osu874GfBtI4K8yc/NhtrEJ2AQwNjb2s7fccstSn8MR++CX5gC48uw1pR5ndnaWl7/85VVEGipzVq8tWcvmrGpf6ee88857YJivPA27hi1n/arKShnDy6ktWc1ZvdI1LDP7XoBTgUe6TH8vcBtFU9dl/knF9fHAQ8AvLGV7p59+eg7Tm6/7Yr75ui+Wfpy77rqrfJhlYM7qtSVr2ZxV7Sv9APfnEmrDoJflrGHDrl9VWSljeDm1Jas5q1e2hg18NmVEXAq8HviNIki3Rm9vcb2/KHhnDbo9SaqSNUxSUwzUjEXEBcDvA7+amc/3WOaYiDh24TbwOuCRbstK0nKyhklqkqV8tcVNwD3ARETsiYi3AdcAxwJ3FKd8X1cse1JEbC1WPQH4QkQ8BNwHfDozPzOUZyFJPVjDJDXdUf0WyMxLukz+WI9l9wIbittPAj9dKp0klWQNk9R0fgO/JElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqUd9mLCKuj4j9EfFIx7RXRsQdEfF4cf2KHuteEBG7ImJ3RFxRZfA6bdk+zblX3clbP/Ndzr3qTrZsn647knQIx+gPWMMO5dhQ063EMbqUV8ZuAC5YNO0K4HOZuQ74XHH/EBGxCvgwcCFwBnBJRJxRKm0DbNk+zZW37mR6Zg6A6Zk5rrx154oYLGoHx+hL3IA1DHBsqPlW6hjt24xl5t3AtxZNvgi4sbh9I7Cxy6pnAbsz88nM/D7wyWK9Vrt62y7mDhw8ZNrcgYNcvW1XTYmkQzlGD2UN+wHHhppupY7RowZc74TM3AeQmfsi4vguy4wD3+i4vwc4u9cDRsQmYBPA2NgYU1NTA0brb6bouAfZxkK33m36MDOXMTs729hsndqSE5qdtcoxWmZfabhKa9hy1q8yrF/D1ZasTc7ZxjFahUGbsaWILtOy18KZuRnYDDAxMZGTk5NDigXX7roHgMnJc4543fF77+w6WMbXrmGYmcuYmppqbLZObckJzc5a5Rgts6+MgCXXsOWsX2VYv4arLVmbnLONY7QKg55N+XREnAhQXO/vsswe4JSO+ycDewfcXmNcfv4Ea1avOmTamtWruPz8iZoSSYdyjC7Jiqxhjg013Uodo4M2Y7cDlxa3LwU+1WWZLwPrIuK0iDgauLhYr9U2rh/ng298DUevmv/Rja9dwwff+Bo2rh+vOZk0zzG6JCuyhjk21HQrdYz2fZsyIm4CJoHjImIP8H7gKuCWiHgb8BTwpmLZk4CPZuaGzHwhIt4JbANWAddn5qPDeRrLa+P6cW667ylmZmbY9vu/WHcc6SUcoz9gDTuUY0NNtxLHaN9mLDMv6THrl7osuxfY0HF/K7B14HSSVJI1TFLT+Q38kiRJNbIZkyRJqpHNmCRJUo1sxiRJkmpkMyZJklQjmzFJkqQa2YxJkiTVyGZMkiSpRjZjkiRJNbIZkyRJqpHNmCRJUo1sxiRJkmpkMyZJklQjmzFJkqQa2YxJkiTVaOBmLCImImJHx+U7EXHZomUmI+LZjmX+sHRiSZKkEXLUoCtm5i7gtQARsQqYBm7rsujnM/P1g25HkoYhIiaAmzsmvQr4w8z8845lJoFPAV8rJt2amX+8TBElrRADN2OL/BLwRGb+c0WPJ0lD5QGlpKao6jNjFwM39Zh3TkQ8FBH/GBE/WdH2JKlKHlBKqk3pV8Yi4mjgV4Eru8x+EPjxzJyNiA3AFmBdj8fZBGwCGBsbY2pqqmy0nmZm5gBKbWNmZo6DBw8ONWdVZmdnzVmxNmStYoxWsa+0RN8DSmAv8J7MfHTxAstZv6pg/RqOtmRtQ842jdEqVPE25YXAg5n59OIZmfmdjttbI+IvI+K4zHymy7Kbgc0AExMTOTk5WUG07q7ddQ8Ak5PnlHqMmZkZhpmzKlNTU+asWBuyVjFGq9hXmq6KA8rlrF9VsH4NR1uytiFnm8ZoFap4m/ISehxRRsSPRkQUt88qtvfNCrYpSVU57AFlZs4Wt7cCqyPiuOUOKGm0lXplLCJeBvx74O0d094BkJnXAb8G/HZEvADMARdnZpbZpiRV7LAHlMDTmZkeUEoallLNWGY+D/zIomnXddy+BrimzDYkaVg8oJTUBFV9tYUktY4HlJKawH+HJEmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjUo1YxHx9YjYGRE7IuL+LvMjIv4iInZHxMMR8TNltidJVbKGSWqCoyp4jPMy85ke8y4E1hWXs4Fri2tJagprmKRaDfttyouAj+e8e4G1EXHikLcpSVWxhkkaurKvjCXw2YhI4K8yc/Oi+ePANzru7ymm7Vv8QBGxCdgEMDY2xtTUVMlovc3MzAGU2sbMzBwHDx4cas6qzM7OmrNibchaxRitYl9puEpq2HLWrypYv4ajLVnbkLNNY7QKZZuxczNzb0QcD9wREV/NzLs75keXdbLbAxVFcDPAxMRETk5OlozW27W77gFgcvKcUo8xMzPDMHNWZWpqypwVa0PWKsZoFftKw1VSw5azflXB+jUcbcnahpxtGqNVKPU2ZWbuLa73A7cBZy1aZA9wSsf9k4G9ZbYpSVWxhklqgoGbsYg4JiKOXbgNvA54ZNFitwNvKc5I+jng2cx8yVuUkrTcrGGSmqLM25QnALdFxMLj/E1mfiYi3gGQmdcBW4ENwG7geeA/lYsrSZWxhklqhIGbscx8EvjpLtOv67idwO8Mug1JGhZrmKSm8Bv4JUmSamQzJkmSVCObMUmSpBrZjEmSJNXIZkySJKlGNmOSJEk1shmTJEmqkc2YJElSjWzGJEmSamQzJkmSVCObMUmSpBqV+Ufh0kC2bJ/m6m272Dszx0lr13D5+RNsXD9+xMtIUh361Sfrl46UzZiW1Zbt01x5607mDhwEYHpmjitv3QnA2iUsY0GTVKd+Ncz6pUH4NqWW1dXbdv1LkVowd+AgV2/bdUTLSFId+tUn65cGYTOmZbV3Zq7v9KUsI0l16FefrF8ahM2YltVJa9f0nb6UZSSpDv3qk/VLgxi4GYuIUyLiroh4LCIejYjf7bLMZEQ8GxE7issflourtrv8/AnWrF51yLQ1q1dx+fkTR7SMVJY1TIPoV5+sXxpEmQ/wvwC8OzMfjIhjgQci4o7M/Mqi5T6fma8vsR2NkIUPsP7e3z7M9w++yHjHmUZTU4/3XUaqkDVMR6xfDbN+aRADN2OZuQ/YV9x+LiIeA8aBxYVMOsTG9ePcdN9TANz89nMGXkYqwxqmQfWrT9YvHalKvtoiIk4F1gNf6jL7nIh4CNgLvCczH+3xGJuATQBjY2NMTU1VEa2rmeKDlGW2MTMzx8GDB4easyqzs7ONy9ntd7A4ZxW/p2Fp4s90sSrGaJN/B1UqW8OWs35VwfpVXr8a1uR9p6k/005tGqNVKN2MRcTLgb8DLsvM7yya/SDw45k5GxEbgC3Aum6Pk5mbgc0AExMTOTk5WTZaT9fuugeAycnBj1iu3XUPMzMzDDNnVaamphqXs9vvYHHOKn5Pw9LEn+liVYzRJv8OqlJFDVvO+lUF61d5/WpYk/edpv5MO7VpjFah1NmUEbGa+SL2icy8dfH8zPxOZs4Wt7cCqyPiuDLblKSqWMMkNUGZsykD+BjwWGb+WY9lfrRYjog4q9jeNwfdpiRVxRomqSnKvE15LvCbwM6I2FFM+wPgxwAy8zrg14DfjogXgDng4szMEtuUpKpYwyQ1QpmzKb8ARJ9lrgGuGXQbkjQs1jBJTeE38EuSJNWokq+20Ojasn2aq7ftYu/MHCe1/MsLR+m5SOpvlPb5UXoueimbMfW0Zfs0V966k7kDBwGYnpnjylt3ArSuCIzSc5HU3yjt86P0XNSdb1Oqp6u37fqXnX/B3IGDXL1tV02JBjdKz0VSf6O0z4/Sc1F3NmPqaW/xDdJLnd5ko/RcJPU3Svv8KD0XdWczpp5OWrvmiKY32Sg9F0n9jdI+P0rPRd3ZjKmny8+fYM3qVYdMW7N6FZefP1FTosGN0nOR1N8o7fOj9FzUnR/gV08LHwz9vb99mO8ffJHxFp/BM0rPRVJ/o7TPj9JzUXcrqhlbODV4unif/dQrPu2g7mPj+nFuuu8pAG5+e/P+4e2RGKXnMmzuKxoFo7TPj9Jz0UutmGZs8anBCzxFWDqU+4okLa8V85mxbqcGL/AUYekH3FckaXmtmGas3ynAniIszXNfkaTltWKasbUvW33Y+f9mzeHnL9iyfZpzr7qTL33tWzzx7Its2T5dRTypcs889z2eePZFTrvi05x71Z1LHqtV7SuSpKVZEc3Y+7bs5NvPHzjsMt/5fwf6/rFa+CzNwoeaX3gRLrt5B+v/+LM2ZWqMLdunee0ffZYnnvkuL7wIyQ8+79VvnFa1r0iSlm7km7H3bdnJX9/7VN/lXkx41807DvtH5o/+/tGun6X59vMHeNfNO3jflp2lskplvW/LTt518w5m5l7aUM0dOMgHbn+057q/8ZF7KttXJElL18qzKedfoXqYuQMvVvq4yfwrXZfdvGOgdf/63qde8sfsmKNX8V/f8BrPPlNltmyf5gO3P9q14epnZu4Ap17x6dIZjmRfecXLVvP+X/lJ9wFJ6iEyc/CVIy4A/iewCvhoZl61aH4U8zcAzwNvzcwH+z3uD5+4Lk+89M8HziWpXfbdeBnf2/d4LPd2h1HDrF/SylO2hg38NmVErAI+DFwInAFcEhFnLFrsQmBdcdkEXDvo9iSpStYwSU1R5m3Ks4DdmfkkQER8ErgI+ErHMhcBH8/5l9/ujYi1EXFiZu473AOf/Nx+/tvn/7JENElt8pZ6NjuUGmb9klaesjWszAf4x4FvdNzfU0w70mUAiIhNEXF/RNxfItOSvWz1sr8jIrXOy1YHa//VyO4rldWw5a5fkkZLmVfGulXoxR9AW8oy8xMzNwObYf4zF7//8/+lRLTeFj5Q/yvrx4d2IoA0Cs79t6/kE/95/n/gDX1fefKy4Tzu4VVWw5arfklqqJI1rEwztgc4peP+ycDeAZZZFt3O6Nq4fpyNRVN2JGenLTR0X3nsK9zyeA50Vps0TD8U8Otn/xgve/7/HvEYrXJfabhW1TBJo6tMM/ZlYF1EnAZMAxcDv75omduBdxafxTgbeLbf58WORFWnzC/8oTlSa599nD/49cmu80bsj5Yaqt8+MDX1zZ5jdBBHsq+0YB+ovYZJEgCZOfCF+dO9/wl4AnhvMe0dwDuK28H82UpPADuBM5fyuKeffnq2wV133VV3hCUxZ/XakrUtOYH7s0QtGvQyjBpm/apWW3JmtierOatXtoaV+tLXzNwKbF007bqO2wn8TpltSNKwWMMkNcHI/zskSZKkJrMZkyRJqpHNmCRJUo1sxiRJkmpU6h+FD0tEPAfsqjvHEhwHPFN3iCUwZ/XakrUtOScy89i6Q1TB+lW5tuSE9mQ1Z/VK1bBSZ1MO0a7MPLPuEP1ExP3mrE5bckJ7srYpZ90ZKmT9qlBbckJ7spqzemVrmG9TSpIk1chmTJIkqUZNbcY21x1gicxZrbbkhPZkNefya8tzMWf12pLVnNUrlbWRH+CXJElaKZr6ypgkSdKKYDMmSZJUo1qbsYh4U0Q8GhEvRsSZHdNPjYi5iNhRXK7rmPezEbEzInZHxF9ERNSVs5h3ZZFlV0ScX2fOLrk/EBHTHT/HDf1y1yUiLiiy7I6IK+rO0ykivl78LncsnL4cEa+MiDsi4vHi+hU1Zbs+IvZHxCMd03pmq+v33iNna8ZnN22pX4fLWsxrZA1r0/iwfg2czfq1IDNruwA/AUwAU8CZHdNPBR7psc59wDlAAP8IXFhjzjOAh4AfBk4DngBW1ZWzS+4PAO/pMr1n7prGwaoiw6uAo4tsZ9Q5Nhfl+zpw3KJpfwpcUdy+AvjvNWX7BeBnOveXXtnq/L33yNmK8XmY59SK+tUna2NrWFvGh/WrVDbrV3Gp9ZWxzHwsM5f8TdURcSLwrzPznpx/1h8HNg4r34LD5LwI+GRmfi8zvwbsBs6qK+cR6Jq7xjxnAbsz88nM/D7wySJjk10E3FjcvpGafr+ZeTfwrUWTe2Wr7ffeI2cvTRufXbWlfsHI1bCmjQ/r14CsXz/Q5M+MnRYR2yPi/0TEzxfTxoE9HcvsKabVZRz4Rsf9hTxNyvnOiHi4eJl14eXeXrnr0rQ8iyXw2Yh4ICI2FdNOyMx9AMX18bWle6le2Zr4c27D+BxEG+oXNL+GtWF8NC3PYtav4alsfA793yFFxP8GfrTLrPdm5qd6rLYP+LHM/GZE/CywJSJ+kvmXyxer5Ls5BszZK8/Qcr4kwGFyA9cCf1Js+0+ADwG/tZz5lqhpeRY7NzP3RsTxwB0R8dW6Aw2oaT/nxo/PttQvaGcNs34tC+vXcFQ6PofejGXmLw+wzveA7xW3H4iIJ4DTme8wT+5Y9GRgb105izyndMkztJyLLTV3RHwE+Ifibq/cdWlankNk5t7ien9E3Mb8S85PR8SJmbmveEtnf60hD9UrW6N+zpn59MLtpo7PttSvYlutq2HWr+Gzfg1H1fWrkW9TRsRYRKwqbr8KWAc8Wbxk+VxE/FxxZs9bgF5HfMvhduDiiPjhiDityHlfU3IWA3nBG4CFM0G65l7ufB2+DKyLiNMi4mjg4iJj7SLimIg4duE28Drmf463A5cWi11KveNwsV7ZGvV7b9H4PCItql/Q4BrWovFh/arWyqxfy3EmQq9L8QT2MH8U+TSwrZj+H4BHmT8j4UHgVzrWObN40k8A11D8F4E6chbz3ltk2UXH2UZ15OyS+38BO4GHiwFyYr/cNY6FDcA/FZneW3eejlyvKsbhQ8WYfG8x/UeAzwGPF9evrCnfTcy/LXagGKNvO1y2un7vPXK2Znz2eE6tqF+Hy3q4n3XdNaxN48P6NXA+61dx8d8hSZIk1aiRb1NKkiStFDZjkiRJNbIZkyRJqpHNmCRJUo1sxiRJkmpkMyZJklQjmzFJkqQa/X+7+gYNl9JBugAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "\n",
    "# 生成正弦波\n",
    "def getSin(amp,freq,phase,sampleList):\n",
    "    return amp*np.sin(-2*math.pi*freq*sampleList+phase)\n",
    "\n",
    "# 生成余弦波\n",
    "def getCos(amp,freq,phase,sampleList):\n",
    "    return amp*np.cos(-2*math.pi*freq*sampleList+phase)\n",
    "\n",
    "# 去掉低振幅噪音\n",
    "def denoise(arr, thresh):\n",
    "    # arr[arr < thresh] = 0\n",
    "    # return arr  \n",
    "    mask= (arr > thresh)\n",
    "    mask= mask + 0 # 用mask加0把True和False转化成1和0显示\n",
    "    # 掩模的值是1，相乘保留数值。反之是0，相乘则为0\n",
    "    return mask * arr\n",
    "\n",
    "# 1. 获得混合波形\n",
    "srate=3000\n",
    "t=np.linspace(0,1,srate)\n",
    "s1=getSin(amp=1.5, freq=30, phase=0,sampleList=t)    \n",
    "s2=getCos(amp=3, freq=5, phase=0,sampleList=t)  \n",
    "s3=getSin(amp=10, freq=100, phase=0,sampleList=t) \n",
    "s4=getCos(amp=20, freq=120, phase=0,sampleList=t) \n",
    "\n",
    "m=s1+s2+s3+s4\n",
    "\n",
    "# #####\n",
    "# # 2. 获得绕线的傅里叶系数\n",
    "# fCoefs=np.zeros(srate,dtype='complex') \n",
    "# # 循环所有可能的频率\n",
    "# for f in range(srate):\n",
    "#     p=m*np.exp(-1j*2*math.pi*f*t) # 让飞轮旋转时间t\n",
    "#     fCoefs[f]=np.sum(p)      # 记录傅里叶系数  \n",
    "# #####\n",
    "\n",
    "# 2.获得傅里叶系数\n",
    "# 第一个参数是混合波形，第二个参数是“飞轮”测试频率上限（即采样率）\n",
    "fCoefs = np.fft.fft(m,srate)\n",
    "\n",
    "# 3. 获得振幅列表：每一个绕线的重心到原点的距离\n",
    "amp_list=2 * np.abs(fCoefs/srate)\n",
    "\n",
    "# 把频率轴从0~300 转变成 0~149 然后 -150~-1\n",
    "freqs = np.fft.fftfreq(len(amp_list), 1/srate)\n",
    "\n",
    "# 然后把 频率轴 和 数据 都变成 0hz 在中间，向左是负频率，向右是正频率的形式\n",
    "amp_shifted=np.fft.fftshift(amp_list)\n",
    "freq_shift=np.fft.fftshift(freqs)\n",
    "\n",
    "#去噪前\n",
    "fg ,ax = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "ax[0].stem(freq_shift,amp_shifted)\n",
    "ax[0].set_xlim([-150,150])\n",
    "ax[0].grid()  # 添加网格线\n",
    "\n",
    "\n",
    "#去噪\n",
    "amp_shifted[(freq_shift>110) | (freq_shift<-110)]= 0  #去除频率大于110hz和小于-110hz的噪音\n",
    "amp_shifted = denoise(amp_shifted, 1) #将振幅小于1的噪声全部去除   \n",
    "\n",
    "#去噪后\n",
    "ax[1].stem(freq_shift,amp_shifted)\n",
    "ax[1].set_xlim([-150,150])\n",
    "ax[1].set_ylim([-1,21])\n",
    "ax[1].grid()  # 添加网格线\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "e18aeb5c-674e-47bc-baa3-440be03a147f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAmMAAAD4CAYAAACg9uHUAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAcO0lEQVR4nO3dfYxldZ3n8fd3G3B7gU37UCAUzICxqQyrq4wEJURTjrMCHXdg3NGFbEbcMWl1NRkTZYQhmTExWdklzu5OcBt7lICzs4CZgZYde21ZoVaNKPJow2BpyzjaD6GjTimltQjtd//oU1pd3Ft1+55z63fOrfcrual7f+fpc6p+58uX+9SRmUiSJKmMf1I6gCRJ0npmMyZJklSQzZgkSVJBNmOSJEkF2YxJkiQVdEzpAL1s2rQpX/rSl5aOsaqf/vSnHH/88aVjrMqczetK1q7kfOCBB36QmROlczTB+tWsruSE7mQ1Z/Pq1rBWNmMnn3wy999/f+kYq5qZmWF6erp0jFWZs3ldydqVnBHxD6UzNMX61ayu5ITuZDVn8+rWMF+mlCRJKshmTJIkqSCbMUmSpIJsxiRJkgqyGZMkSSpo1WYsIk6PiHsi4vGIeCwi/rAaf0FE3BUR365+Pr/P9hdFxGxE7ImIq5o+gaOx46F9XHDt3Zx51We54Nq72fHQvpJxpFYbh+tlnOqXpPE1yDNjzwLvz8zfAF4DvCcizgauAr6QmZuBL1SPjxARG4CPARcDZwOXV9uuuR0P7ePq23ezb26BBPbNLXD17bs7+R8YadTG6HoZi/olabyt2oxl5oHMfLC6/xTwODAJXALcXK12M3Bpj83PA/Zk5hOZ+XPg1mq7NXfdrlkWnjl0xNjCM4e4btdsiThSq43L9TIu9UvSeDuqL32NiDOAc4CvASdn5gE4XPAi4qQem0wC31/yeC/w6j773gpsBZiYmGBmZuZooq1q39xC3/FhjzU/P994zlEwZ/O6knXYnKO4Xkrrcv0ahXGfwyV0Jas522fgZiwiTgD+BnhfZv4kIgbarMdY9loxM7cD2wGmpqay6W/dnfzq3T3/AzO5aePQ3/DblW8HNmfzupJ12JyjuF5K6nr9GoVxn8MldCWrOdtnoE9TRsSxHC5kf5WZt1fDT0bEKdXyU4CDPTbdC5y+5PFpwP7h4w7vygun2HjshiPGNh67gSsvnCoRR2q1cbpexqF+SRpvg3yaMoBPAo9n5p8tWXQncEV1/wrgMz02/zqwOSLOjIjjgMuq7dbcpedM8pE3v5zjNhw+5clNG/nIm1/OpedMlogjtdq4XC/jUr8kjbdBXqa8APh9YHdEPFyN/TFwLfDpiHgH8D3gLQARcSrwiczckpnPRsR7gV3ABuDGzHys4XMY2KXnTHLLfd8D4LZ3nl8qhtQJY3K9jE39kjS+Vm3GMvPL9H7vBMAbeqy/H9iy5PFOYOewASVpWNYvSV3gN/BLkiQVZDMmSZJUkM2YJElSQTZjkiRJBdmMSZIkFWQzJkmSVJDNmCRJUkE2Y5IkSQXZjEmSJBVkMyZJklSQzZgkSVJBNmOSJEkF2YxJkiQVZDMmSZJUkM2YJElSQTZjkiRJBR2z2goRcSPwJuBgZr6sGrsNmKpW2QTMZeYre2z7XeAp4BDwbGae20hqSRqQNUxS263ajAE3AdcDn1ocyMx/u3g/Ij4K/HiF7V+fmT8YNqAk1XQT1jBJLbZqM5aZX4yIM3oti4gA3gr8VsO5JKkR1jBJbTfIM2MreS3wZGZ+u8/yBD4fEQl8PDO399tRRGwFtgJMTEwwMzNTM1pvc3MLAI3sf35+fmQ5m2TO5nUla92cTV4vLdVIDVur+tWk9TKH11JXspqzhTJz1RtwBvBoj/FtwPtX2O7U6udJwCPA6wY53llnnZWj8tYbvpJvveErjezrnnvuaWQ/o2bO5nUla92cTV4vKwHuzwFqw7C3taxho6xfTVovc3gtdSWrOZtXt4YN/WnKiDgGeDNw2wqN3v7q50HgDuC8YY8nSU2yhklqizpfbfHbwDczc2+vhRFxfEScuHgfeCPwaI3jSVKTrGGSWmHVZiwibgHuBaYiYm9EvKNadBlwy7J1T42IndXDk4EvR8QjwH3AZzPzc81Fl6TVWcMktd0gn6a8vM/423uM7Qe2VPefAF5RM58k1WINk9R2fgO/JElSQTZjkiRJBdmMSZIkFWQzJkmSVJDNmCRJUkE2Y5IkSQXZjEmSJBVkMyZJklSQzZgkSVJBNmOSJEkF2YxJkiQVZDMmSZJUkM2YJElSQTZjkiRJBdmMSZIkFWQzJkmSVNCqzVhE3BgRByPi0SVjH4qIfRHxcHXb0mfbiyJiNiL2RMRVTQaXpEFYwyS13SDPjN0EXNRj/L9k5iur287lCyNiA/Ax4GLgbODyiDi7TlhJGsJNWMMktdiqzVhmfhH40RD7Pg/Yk5lPZObPgVuBS4bYjyQNzRomqe2OqbHteyPibcD9wPsz8x+XLZ8Evr/k8V7g1f12FhFbga0AExMTzMzM1IjW39zcAkAj+5+fnx9ZziaZs3ldyVo3Z5PXSws1VsPWqn41ab3M4bXUlazmbJ9hm7FtwIeBrH5+FPiDZetEj+2y3w4zczuwHWBqaiqnp6eHjLaybbP3AjA9fX7tfc3MzDCqnE0yZ/O6krVuziavl5ZptIatVf1q0nqZw2upK1nN2T5DfZoyM5/MzEOZ+QvgLzj8dP5ye4HTlzw+Ddg/zPEkqUnWMEltMlQzFhGnLHn4u8CjPVb7OrA5Is6MiOOAy4A7hzmeJDXJGiapTVZ9mTIibgGmgRdFxF7gT4HpiHglh5+y/y7wzmrdU4FPZOaWzHw2It4L7AI2ADdm5mOjOAlJ6scaJqntVm3GMvPyHsOf7LPufmDLksc7ged8ZFyS1oo1TFLb+Q38kiRJBdmMSZIkFWQzJkmSVJDNmCRJUkE2Y5IkSQXZjEmSJBVkMyZJklSQzZgkSVJBNmOSJEkF2YxJkiQVZDMmSZJUkM2YJElSQTZjkiRJBdmMSZIkFWQzJkmSVNCqzVhE3BgRByPi0SVj10XENyPiGxFxR0Rs6rPtdyNid0Q8HBH3N5hbkgZiDZPUdoM8M3YTcNGysbuAl2XmvwS+BVy9wvavz8xXZua5w0WUpFpuwhomqcVWbcYy84vAj5aNfT4zn60efhU4bQTZJKk2a5iktovMXH2liDOAv83Ml/VY9r+A2zLzf/RY9vfAPwIJfDwzt69wjK3AVoCJiYlXffrTnx70HI7KR762AMDVr95Ye1/z8/OccMIJtfczauZsXley1s3Z5PWykte//vUPjPKZp1HXsLWqX01aL3N4LXUlqzmbV7uGZeaqN+AM4NEe49cAd1A1dT2Wn1r9PAl4BHjdIMc766yzclTeesNX8q03fKWRfd1zzz2N7GfUzNm8rmStm7PJ62UlwP05QG0Y9raWNWyU9atJ62UOr6WuZDVn8+rWsKE/TRkRVwBvAv5dFaRXo7e/+nmwKnjnDXs8SWqSNUxSWwzVjEXERcAHgd/JzJ/1Wef4iDhx8T7wRuDRXutK0lqyhklqk0G+2uIW4F5gKiL2RsQ7gOuBE4G7qo9831Cte2pE7Kw2PRn4ckQ8AtwHfDYzPzeSs5CkPqxhktrumNVWyMzLewx/ss+6+4Et1f0ngFfUSidJNVnDJLWd38AvSZJUkM2YJElSQTZjkiRJBdmMSZIkFWQzJkmSVJDNmCRJUkE2Y5IkSQXZjEmSJBVkMyZJklSQzZgkSVJBNmOSJEkF2YxJkiQVZDMmSZJUkM2YJElSQTZjkiRJBdmMSZIkFbRqMxYRN0bEwYh4dMnYCyLiroj4dvXz+X22vSgiZiNiT0Rc1WTwUnY8tI8Lrr2bM6/6LO+f+Rk7HtpXOpL0HIvz9O2f+ykXXHv3up6n1rBfsX6pC9Zj/RrkmbGbgIuWjV0FfCEzNwNfqB4fISI2AB8DLgbOBi6PiLNrpS1sx0P7uPr23eybWyCBH/6/5Orbd6+LiaLuWDpPAfbNLaz3eXoT1jDrlzphvdavVZuxzPwi8KNlw5cAN1f3bwYu7bHpecCezHwiM38O3Fpt11nX7Zpl4ZlDR4wtPHOI63bNFkokPZfz9EjWsMOcF+qC9TpPjxlyu5Mz8wBAZh6IiJN6rDMJfH/J473Aq/vtMCK2AlsBJiYmmJmZGTLayuaqbnuY/S926r3GR5W3CfPz863Ot6grOaHdWZucp3Wul5ZrtIatVf2qw/o1el3J2uacXZ2ndQ3bjA0ieoxlv5UzczuwHWBqaiqnp6dHEmrb7L0ATE+ff9TbTn717p4TZXLTRkaVtwkzMzOtzreoKzmh3VmbnKd1rpcxMHANW6v6VYf1a/S6krXNObs6T+sa9tOUT0bEKQDVz4M91tkLnL7k8WnA/iGP1wpXXjjFxmM3HDG28dgNXHnhVKFE0nM5Twey7mqY80JdsF7n6bDN2J3AFdX9K4DP9Fjn68DmiDgzIo4DLqu266xLz5nkI29+OcdtOPxre+E/DT7y5pdz6TmThZNJv7J8nk5u2ug8fa51V8OsX+qC9Vq/Vn2ZMiJuAaaBF0XEXuBPgWuBT0fEO4DvAW+p1j0V+ERmbsnMZyPivcAuYANwY2Y+NprTWDuXnjPJLfd9D4B3Tz3N9JhPEHXT4jydm5tj1wd/q3Scoqxhv2L9Uhesx/q1ajOWmZf3WfSGHuvuB7YsebwT2Dl0OkmqyRomqe38Bn5JkqSCbMYkSZIKshmTJEkqyGZMkiSpIJsxSZKkgmzGJEmSCrIZkyRJKshmTJIkqSCbMUmSpIJsxiRJkgqyGZMkSSrIZkySJKkgmzFJkqSCbMYkSZIKshmTJEkqyGZMkiSpoKGbsYiYioiHl9x+EhHvW7bOdET8eMk6f1I7sSQ1wBomqS2OGXbDzJwFXgkQERuAfcAdPVb9Uma+adjjSNIoWMMktUVTL1O+AfhOZv5DQ/uTpLVkDZNUzNDPjC1zGXBLn2XnR8QjwH7gA5n5WK+VImIrsBVgYmKCmZmZhqIdaW5uAaDW/hf3MT9/aGQ5mzQ/P2/OhnUh69zcAocO1ZujTVwvHVGrhq1V/WqC9Wt0upK1CzmbqF9dUrsZi4jjgN8Bru6x+EHg1zNzPiK2ADuAzb32k5nbge0AU1NTOT09XTdaT9tm7wVgevr82vs44YSnGVXOJs3MzJizYV3Ium32Xubm5mrlbOJ6absmatha1a8mWL9GpytZu5CzifrVJU28THkx8GBmPrl8QWb+JDPnq/s7gWMj4kUNHFOSmmINk1RUE83Y5fR5ej8iXhwRUd0/rzreDxs4piQ1xRomqahaL1NGxD8D/hXwziVj7wLIzBuA3wPeHRHPAgvAZZmZdY4pSU2xhklqg1rNWGb+DHjhsrEblty/Hri+zjEkaVSsYZLawG/glyRJKshmTJIkqSCbMUmSpIJsxiRJkgqyGZMkSSrIZkySJKkgmzFJkqSCbMYkSZIKshmTJEkqyGZMkiSpIJsxSZKkgmzGJEmSCrIZkyRJKshmTJIkqSCbMUmSpIJqNWMR8d2I2B0RD0fE/T2WR0T8eUTsiYhvRMRv1jmeJDXJGiapDY5pYB+vz8wf9Fl2MbC5ur0a2Fb9lKS2sIZJKmrUL1NeAnwqD/sqsCkiThnxMSWpKdYwSSNX95mxBD4fEQl8PDO3L1s+CXx/yeO91diBmseVpCY0UsMiYiuwFWBiYoKZmZmRBa5rbm4BgPn5Q63OuWh+fr4TOaE7WbuQc25ugUOHujFHm1C3GbsgM/dHxEnAXRHxzcz84pLl0WOb7LWjtSpmi4Wozv4tZqPRlZzQjaxNFLMmrpeWa6SGVU3cdoCpqamcnp4eSdgmbJu9F4ATTniaNudcNDMz04mc0J2sXci5bfZe5ubmWp+zKbWasczcX/08GBF3AOcBSwvZXuD0JY9PA/b32deaFLPFQjQ9fX7tfVjMmtWVnNCNrE0UsyaulzZrsoZJ0rCGfs9YRBwfEScu3gfeCDy6bLU7gbdVn0h6DfDjzPQlSknFWcMktUWdZ8ZOBu6IiMX9/M/M/FxEvAsgM28AdgJbgD3Az4B/Xy+uJDXGGiapFYZuxjLzCeAVPcZvWHI/gfcMewxJGhVrmKS28Bv4JUmSCrIZkyRJKshmTJIkqSCbMUmSpIJsxiRJkgqyGZMkSSrIZkySJKkgmzFJkqSCbMYkSZIKshmTJEkqyGZMkiSpoDr/ULh01HY8tI/rds2yf26BUzdt5MoLp7j0nMmh15OktTRIbbJ+6WjZjGnN7HhoH1ffvpuFZw4BsG9ugatv3w1wRKEadD1JWkuD1Cbrl4bhy5RaM9ftmv1lgVq08Mwhrts1O9R6krSWBqlN1i8Nw2ZMa2b/3MJA44OuJ0lraZDaZP3SMGzGtGZO3bRxoPFB15OktTRIbbJ+aRhDN2MRcXpE3BMRj0fEYxHxhz3WmY6IH0fEw9XtT+rFVZddeeEUG4/dcMTYxmM3cOWFU0OtJ9VhDdPRGqQ2Wb80jDpv4H8WeH9mPhgRJwIPRMRdmfl3y9b7Uma+qcZxNCYW37z6R3/9DX5+6BdM9vmU0aDrSTVZw3RUBqlN1i8NY+hmLDMPAAeq+09FxOPAJLC8kEm/dOk5k9xy3/cAuO2d59deTxqWNUzDGKQ2Wb90tBr5aouIOAM4B/haj8XnR8QjwH7gA5n5WJ99bAW2AkxMTDAzM9NEtOeYq95EWWf/i/uYnz80spxNmp+fb1XOfn+D5Tmb+FuNStt+p73MzS1w6FC9Odrmv0GT6tawtapfTbB+1TdIDWvztdPG3+lyTdSvLqndjEXECcDfAO/LzJ8sW/wg8OuZOR8RW4AdwOZe+8nM7cB2gKmpqZyenq4bradts/cCMD09/P+tLO7jhBOeZlQ5mzQzM9OqnP3+BstzNvG3GpW2/U572TZ7L3Nzc7Vytvlv0JQmatha1a8mWL/qG6SGtfnaaePvdLkm6leX1Po0ZUQcy+Ei9leZefvy5Zn5k8ycr+7vBI6NiBfVOaYkNcUaJqkN6nyaMoBPAo9n5p/1WefF1XpExHnV8X447DElqSnWMEltUedlyguA3wd2R8TD1dgfA78GkJk3AL8HvDsingUWgMsyM2scU5KaYg2T1Ap1Pk35ZSBWWed64PphjyFJo2INk9QWfgO/JElSQY18tYXG046H9nHdrln2zy1w6hh8ceG4nY+klY3TNT9O56LnshlTTzse2sfVt+9m4ZlDAOybW+Dq23cDdLIAjNv5SFrZOF3z43Qu6s2XKdXTdbtmf3nhL1p45hDX7ZotlKiecTsfSSsbp2t+nM5FvdmMqaf91bdHDzreduN2PpJWNk7X/Didi3qzGVNPp27aeFTjbTdu5yNpZeN0zY/Tuag3mzH1dOWFU2w8dsMRYxuP3cCVF04VSlTPuJ2PpJWN0zU/Tuei3nwDv3pafFPoH/31N/j5oV8w2fFP74zb+Uha2Thd8+N0LurNZkx9XXrOJLfc9z0Abntn+/6x26M1bucjaWXjdM2P07nouXyZUpIkqSCbMUmSpIJsxiRJkgqyGZMkSSrIZkySJKkgmzFJkqSCbMYkSZIKiswcfuOIi4D/BmwAPpGZ1y5bHtXyLcDPgLdn5oOr7fd5p2zOU674r0PnktQtB25+H08f+Has9XFHUcOsX9L6U7eGDf3MWERsAD4GXAycDVweEWcvW+1iYHN12wpsG/Z4ktQka5iktqjzDfznAXsy8wmAiLgVuAT4uyXrXAJ8Kg8//fbViNgUEadk5oGVdnzaUwf5j1/67zWiSeqSt5U57EhqmPVLWn/q1rA67xmbBL6/5PHeauxo1wEgIrZGxP0RcX+NTJI0qMZqmPVLUh11nhnr9dro8jegDbLO4cHM7cB2OPyeiw++9j/UiCapU554X4mjNlbDrF/SOlezhtV5ZmwvcPqSx6cB+4dYR5JKsIZJaoU6zdjXgc0RcWZEHAdcBty5bJ07gbfFYa8Bfrza+8UkaY1YwyS1Q2YOfePwx72/BXwHuKYaexfwrup+cPjTSt8BdgPnDrLfs846K7vgnnvuKR1hIOZsXleydiUncH/WqEXD3kZRw6xfzepKzszuZDVn8+rWsDrvGSMzdwI7l43dsOR+Au+pcwxJGhVrmKQ28Bv4JUmSCrIZkyRJKshmTJIkqSCbMUmSpIJq/UPhoxIRTwGzpXMM4EXAD0qHGIA5m9eVrF3JOZWZJ5YO0QTrV+O6khO6k9WczatVw2p9mnKEZjPz3NIhVhMR95uzOV3JCd3J2qWcpTM0yPrVoK7khO5kNWfz6tYwX6aUJEkqyGZMkiSpoLY2Y9tLBxiQOZvVlZzQnazmXHtdORdzNq8rWc3ZvFpZW/kGfkmSpPWirc+MSZIkrQs2Y5IkSQUVbcYi4i0R8VhE/CIizl0yfkZELETEw9XthiXLXhURuyNiT0T8eUREqZzVsqurLLMRcWHJnD1yfygi9i35PW5ZLXcpEXFRlWVPRFxVOs9SEfHd6m/58OLHlyPiBRFxV0R8u/r5/ELZboyIgxHx6JKxvtlK/d375OzM/OylK/VrpazVslbWsC7ND+vX0NmsX4sys9gN+A1gCpgBzl0yfgbwaJ9t7gPOBwL438DFBXOeDTwCPA84E/gOsKFUzh65PwR8oMd439yF5sGGKsNLgOOqbGeXnJvL8n0XeNGysf8MXFXdvwr4T4WyvQ74zaXXS79sJf/ufXJ2Yn6ucE6dqF+rZG1tDevK/LB+1cpm/apuRZ8Zy8zHM3Pgb6qOiFOAf56Z9+bhs/4UcOmo8i1aIeclwK2Z+XRm/j2wBzivVM6j0DN3wTznAXsy84nM/Dlwa5WxzS4Bbq7u30yhv29mfhH40bLhftmK/d375OynbfOzp67ULxi7Gta2+WH9GpL161fa/J6xMyPioYj4vxHx2mpsEti7ZJ291Vgpk8D3lzxezNOmnO+NiG9UT7MuPt3bL3cpbcuzXAKfj4gHImJrNXZyZh4AqH6eVCzdc/XL1sbfcxfm5zC6UL+g/TWsC/OjbXmWs36NTmPzc+T/HFJE/B/gxT0WXZOZn+mz2QHg1zLzhxHxKmBHRPwLDj9dvlwj380xZM5+eUaW8zkBVsgNbAM+XB37w8BHgT9Yy3wDalue5S7IzP0RcRJwV0R8s3SgIbXt99z6+dmV+gXdrGHWrzVh/RqNRufnyJuxzPztIbZ5Gni6uv9ARHwHOIvDHeZpS1Y9DdhfKmeV5/QeeUaWc7lBc0fEXwB/Wz3sl7uUtuU5Qmbur34ejIg7OPyU85MRcUpmHqhe0jlYNOSR+mVr1e85M59cvN/W+dmV+lUdq3M1zPo1etav0Wi6frXyZcqImIiIDdX9lwCbgSeqpyyfiojXVJ/seRvQ7//41sKdwGUR8byIOLPKeV9bclYTedHvAoufBOmZe63zLfF1YHNEnBkRxwGXVRmLi4jjI+LExfvAGzn8e7wTuKJa7QrKzsPl+mVr1d+9Q/PzqHSofkGLa1iH5of1q1nrs36txScR+t2qE9jL4f+LfBLYVY3/G+AxDn8i4UHgXy/Z5tzqpL8DXE/1rwiUyFktu6bKMsuSTxuVyNkj918Cu4FvVBPklNVyF5wLW4BvVZmuKZ1nSa6XVPPwkWpOXlONvxD4AvDt6ucLCuW7hcMviz1TzdF3rJSt1N+9T87OzM8+59SJ+rVS1pV+16VrWJfmh/Vr6HzWr+rmP4ckSZJUUCtfppQkSVovbMYkSZIKshmTJEkqyGZMkiSpIJsxSZKkgmzGJEmSCrIZkyRJKuj/A627ukzWfMFrAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 720x288 with 2 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import math\n",
    "\n",
    "# 生成正弦波\n",
    "def getSin(amp,freq,phase,sampleList):\n",
    "    return amp*np.sin(-2*math.pi*freq*sampleList+phase)\n",
    "\n",
    "# 生成余弦波\n",
    "def getCos(amp,freq,phase,sampleList):\n",
    "    return amp*np.cos(-2*math.pi*freq*sampleList+phase)\n",
    "\n",
    "# 去掉低振幅噪音\n",
    "def denoise(arr, thresh):\n",
    "    # arr[arr < thresh] = 0\n",
    "    # return arr  \n",
    "    mask= (arr > thresh)\n",
    "    mask= mask + 0 # 用mask加0把True和False转化成1和0显示\n",
    "    # 掩模的值是1，相乘保留数值。反之是0，相乘则为0\n",
    "    return mask * arr\n",
    "\n",
    "# 1. 获得混合波形\n",
    "srate=3000\n",
    "t=np.linspace(0,1,srate)\n",
    "s1=getSin(amp=1.5, freq=30, phase=0,sampleList=t)    \n",
    "s2=getCos(amp=3, freq=5, phase=0,sampleList=t)  \n",
    "s3=getSin(amp=10, freq=100, phase=0,sampleList=t) \n",
    "s4=getCos(amp=20, freq=120, phase=0,sampleList=t) \n",
    "\n",
    "m=s1+s2+s3+s4\n",
    "\n",
    "# 2. 获得绕线的傅里叶系数\n",
    "fCoefs=np.zeros(srate,dtype='complex') \n",
    "# 循环所有可能的频率\n",
    "for f in range(srate):\n",
    "    p=m*np.exp(-1j*2*math.pi*f*t) # 让飞轮旋转时间t\n",
    "    fCoefs[f]=np.sum(p)      # 记录傅里叶系数  \n",
    "\n",
    "# 3. 获得振幅列表：每一个绕线的重心到原点的距离\n",
    "amp_list=2 * np.abs(fCoefs/srate)\n",
    "\n",
    "# 把频率轴从0~300 转变成 0~149 然后 -150~-1\n",
    "freqs = np.fft.fftfreq(len(amp_list), 1/srate)\n",
    "\n",
    "# 然后把 频率轴 和 数据 都变成 0hz 在中间，向左是负频率，向右是正频率的形式\n",
    "amp_shifted=np.fft.fftshift(amp_list)\n",
    "freq_shift=np.fft.fftshift(freqs)\n",
    "\n",
    "#去噪前\n",
    "fg ,ax = plt.subplots(1,2,figsize=(10,4))\n",
    "\n",
    "ax[0].stem(freq_shift,amp_shifted)\n",
    "ax[0].set_xlim([-150,150])\n",
    "ax[0].grid()  # 添加网格线\n",
    "\n",
    "\n",
    "#去噪\n",
    "amp_shifted[(freq_shift>110) | (freq_shift<-110)]= 0  #去除频率大于110hz和小于-110hz的噪音\n",
    "amp_shifted = denoise(amp_shifted, 1) #将振幅小于1的噪声全部去除   \n",
    "\n",
    "#去噪后\n",
    "ax[1].stem(freq_shift,amp_shifted)\n",
    "ax[1].set_xlim([-150,150])\n",
    "ax[1].set_ylim([-1,21])\n",
    "ax[1].grid()  # 添加网格线\n",
    "\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
